Exotic vs. conventional scaling and universality in a 
disordered bilayer quantum Heisenberg antiferromagnet 
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We present Monte-Carlo simulations of a two-dimensional (2d) bilayer quantum Heisenberg an- 
tiferromagnet with random dimer dilution. In contrast to exotic scaling scenarios found in other 
random quantum systems, the quantum phase transition in this system is characterized by a finite- 
disorder fixed point with power-law scaling. After accounting for corrections to scaling, with a 
leading irrelevant exponent of a; « 0.48, we find universal critical exponents z — 1.310(6) and 
1/ = 1.16(3). We discuss the consequences of these findings and suggest new experiments. 
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Quantum phase transitions (QPT) under the influence 
of quenched disorder are a topic of great current interest. 
Experimental examples range from localized and itin- 
erant ^ quantum magnets to heavy-fermion compounds 
01 , high-tcmperature superconductors , and to metal- 
insulator and superconductor- insulator transitions . 
These systems display rich new physics but many are 
still poorly understood. In the context of classical phase 
transitions, the interplay between disorder and critical 
fluctuations has a long history. Harris derived a crite- 
rion for the stability of a critical point against disorder: 
If the correlation length exponent v fulflUs the inequal- 
ity v > 2/d, where d is the spatial dimensionality, the 
critical behavior is not influenced by weak disorder. If 
a clean critical point violates the Harris criterion, the 
generic result of introducing disorder is a new (flnite- 
disorder) critical point with power-law scaling and new 
critical exponents which fulfill the Harris criterion 

At QPTs, order-parameter fluctuations in space and 
time must be considered. Quenched disorder is perfectly 
correlated in time direction. As a result, disorder ef- 
fects at QPTs are generically stronger than at classi- 
cal transitions. Prominent consequences are the infinite- 
randomness critical points in Id random spin chains 9] 
and in Id ^3 and 2d [Oi random quantum Ising 
models. At these critical points, the dynamical scaling is 
activated, i.e., correlation time and correlation length 
f obey In^T- ~ (At conventional critical points, this 
relation is a power law, S,t ~ , with a universal dynami- 
cal exponent z). In itinerant electron systems, the effects 
of impurities can be even more dramatic. For Ising sym- 
metry, the interplay of quenched disorder and Landau 
damping of the order parameter fluctuations completely 
destroys the sharp QPT by smearing . Further exotic 
phenomena include non-universal, continuously varying 
exponents, observed in the Griffiths region associated 
with a QPT [ill [ij or at certain impurity QPTs 
|l5|. On the other hand, the stable low-energy (as op- 
posed to critical) fixed point of random Heisen berg mod- 
els in d > 2 has been shown to be conventional [la| • Pre- 
liminary results for the critical point in these models 



suggest that the infinite-randomness fixed point is unsta- 
ble, but no definite conclusion on the fate of the transition 
has been reached. These results lead to the general ques- 
tion whether all QPTs in presence of quenched disorder 
are unconventional. 

In this Letter, we provide a "proof of principle" that 
this is not the case: The QPT of a dimer-diluted spin-1/2 
bilayer quantum Heisenberg antiferromagnet is shown to 
exhibit a conventional finite-disorder critical point with 
power-law dynamical scaling and universal critical expo- 
nents. After accounting for corrections to scaling char- 
acterized by an irrelevant exponent iv « 0.48 we flnd the 
asymptotic dynamical and correlation length exponents 
to be z = 1.310(6) and i/ = 1.16(3) (fulflUing the Harris 
criterion > 2/d = 1 0,I1)- 

Our starting point is a bilayer quantum Heisenberg an- 
tiferromagnet as depicted in the inset of Fig. The spins 
in each 2d layer interact via nearest neighbor exchange 
J||, and the interplane coupling is J±. The clean version 
of this model has been studied extensively 0, . For 
J± ^ J|| , neighboring spins from the two layers form 
singlets, and the ground state is paramagnetic. In con- 
trast, for J|| ;» J± the system develops Neel order. Both 
phases are separated by a QPT at J±/J\\ ~ 2.525. Ran- 
dom disorder is introduced by removing pairs (dimers) of 
adjacent spins, one from each layer. The Hamiltonian of 
the model with dimer dilution is: 



H 



J 



E 



(1) 



{i.j) 
a=l,2 



and 6^=0 (ei=l) with probability p (1 — p). 

The phase diagram of the dimer-diluted bilayer Heisen- 
berg model has been studied by Sandvik 19] and Vajk 
and Greven |23, see Fig.^ For small J±, magnetic order 
survives up to the percolation threshold Pp « 0.4072, and 
a multicritical point exists at p = Pp and J±/ J\\ ~ 0.16. 
We focus on the generic transition at < p < Pp, driven 
by J±, where the results of Refs. [l^l20j are inconclusive. 

To determine the critical behavior at the QPT effec- 
tively, we proceed by mapping the quantum Hamiltonian 
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FIG. 1: Phase diagram |20ll of the diluted bilayer Heisenberg 
antiferromagnet, as function of J±/J\\ and dilution p. The 
dashed line is the percolation threshold, the open dot is the 
multicritical point of Refs. [iflL l2(il|. The arrow indicates the 
QPT studied here. Inset: The model: Quantum spins (ar- 
rows) reside on the two parallel square lattices. The spins in 
each plane interact with the coupling strength Jy . Interplane 
coupling is J±. Dilution is done by removing dimers. 



onto a classical model. First we note that the low- 
energy properties of bilayer quantum antiferromagnets 
are represented by a (2-|-l)-dimensional 0(3) quantum 
rotor model jl^l with the rotor coordinate n.^ correspond- 
ing to S.i^i — Si^2 and the angular momentum repre- 
senting S.i_i -f Si_2 (see, e.g., chapter 5 of Ref. [l^l)- This 
quantum rotor model in turn is equivalent to a 3d clas- 
sical Heisenberg model with the disorder perfectly corre- 
lated in imaginary time direction, as can be easily seen 
from a path integral representation of the partition func- 
tion. Thus, our classical Hamiltonian reads: 



H = K 



i,T + l J 



(2) 



where r is an 0(3) unit vector. The coupling constant 
PK of the classical model is related to the ratio J\\/J± 
of the quantum model. Here, (3 = 1/T where T is an 
effective "classical" temperature, not equal to the real 
temperature which is zero. We set K = 1 and drive 
the classical system through the transition by tuning the 
classical temperature T. 

As an aside, we note that dimer dilution in the quan- 
tum model does not introduce random Berry phases 
because the Berry phase contributions from the two spins 
of each unit cell cancel 0,1131 • In contrast, for site dilu- 
tion, the physics changes completely: The random Berry 
phases (which have no classical analogue) are equiva- 
lent to impurity- induced moments [j^ l , and those become 
weakly coupled via bulk excitations. Thus, for all p < pp 
the ground state shows long-range order, independent of 
J±/J\\l This effect is absent for dimer dilution, and both 
phases of the clean system survive for small p . 

The classical model is studied by Monte-Carlo 
simulations using the efficient Wolff cluster algorithm 



|25l l26| . We investigate linear sizes up to L = 120 in 
space direction and Lr = 384 in imaginary time, for im- 
purity concentrations p = i , i , | and i . The results are 
averaged over 10'^ - 10^ disorder realizations. Each sam- 
ple is equilibrated using 100 Monte-Carlo sweeps (spin- 
flips per site). For large dilutions, P — j and | we per- 
form both Wolff and Metropolis sweeps to equilibrate 
small dangling clusters. During the measurement period 
of another 100-200 sweeps we calculate magnetization, 
susceptibility, specific heat and correlation functions. 

A quantity particularly suitable to locate the critical 
point and to extract high precision values for the expo- 
nents z and V is the Binder ratio: 

_ r (|M|4) ] 

where M = J2i r ^i-r, [• ■ ^av denotes the disorder aver- 
age and (...) denotes the Monte-Carlo average for each 
sample. This quantity has scale dimension 0. Thus, its 
finite-size scaling form is given by 

gav = ~gc{tL^'\Lr/L'') or (4) 
gav = ~gA{tL^'\\og{Lr)/L^^) (5) 

for conventional scaling or for activated scaling, respec- 
tively. Two important characteristics follow: (i) For fixed 
L, Qav has a peak as a function of Lr- The peak position 
^max ma^j-jjg |;]-^g optimal sample shape, where the ratio 
Lr/L roughly behaves like the corresponding ratio of the 
correlation lengths in time and space directions, ^r/C- 
the critical temperature Tc, the peak value g'^^^ is inde- 
pendent of L. Thus, for power law scaling, plotting gav 
vs. Lr/L"^"-^ at Tc should collapse the data, without the 
need for a value of z. In contrast, for activated scaling 
the gav data should collapse when plotted as a function of 
\og{Lr)/\og{L^^^). (ii) For samples of the optimal shape 
(Lt- = L™°^), plots of gav vs. temperature for different L 
cross at Tc- Based on these two characteristics, we use a 
simple iterative procedure to determine both the optimal 
shapes and the location of the critical point. 

We now turn to our results. To distinguish between 
activated and power-law dynamical scaling we perform 
a series of calculations at the critical temperature. The 
upper panel of Fig.[21shows the Binder ratio gav as a func- 
tion of Lr for various L = 5 . . . 100 and dilution p = ^ 
ai T — Tc — 1.1955. The statistical error of gav is be- 
low 0.1% for the smaller sizes and not more than 0.2% 
for the largest systems. As expected at Tc, the maximum 
Binder ratio for each of the curves does not depend on L- 
To test the conventional power-law scaling form, eq. I@J, 
we plot gavjg^v^ as a function of Lr/L^"-^ in the lower 
panel of Fig. |3 The data scale extremely well, giving 
statistical errors of L^'^^ in the range between 0.3% and 
1%. For comparison, the inset shows a plot of gav as a 
function of log(L^)/ log(L™°'^) corresponding to eq. lO- 
The data clearly do not scale which rules out the acti- 
vated scaling scenario. The results for the other impurity 
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FIG. 2: Upper panel: Binder ratio gav as a function of Lr 
for various L (p = |). Lower panel: Power-law scaling plot 
ffa-u/ff™"" vs. Lt/L™"^ Inset: Activated scaling plot ga^/g™"^ 
vs. j/ = log(LO/log(ir"). 
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FIG. 3: U^°'^ / L vs. L for four disorder concentrations p ~ \, 
s: Fit to 
0.48(3). 



1 I and i. Solid lines: Fit to L^'"' = cii^(l + fcL-"!) with 



=^1.310(6) and 



concentrations p = i, |, i are completely analogous. 

Having established conventional power-law dynamical 
scaling, we proceed to determine the dynamical exponent 
z. In Fig. 131 we plot L™"^ vs. L for all four dilutions p. 
The curves show significant deviations from pure power- 
law behavior which can be attributed to corrections to 
scaling due to irrelevant operators. In such a situation, 
a direct power-law fit of the data will only yield effective 
exponents. To find the true asymptotic exponents we take 
the leading correction to scaling into account by using 
the ansatz = aL^(l + bL-^^) with universal 

(dilution-independent) exponents z and wi but dilution- 
dependent a and h. A combined fit of all four curves 
gives z = 1.310(6) and wi = 0.48(3) where the number in 
brackets is the standard deviation of the last given digit. 
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FIG. 4: Scaling plot of Qav vs. (T - Tc)xl for p = 0.2. is 
the factor necessary to scale the data onto a master curve. 



The fit is of high quality (x^ ~ 0.7) and robust against 
removing complete data sets or removing points form the 
lower or upper end of each set. We thus conclude that 
the asymptotic dynamical exponent z is indeed universal. 
(Note that the leading corrections to scaling vanish very 
close to p = |; the curvature of the L^^^{L) curves in 
Fig. 121 is opposite above and below this concentration.) 

To find the correlation length exponent z/, we perform 
simulations in the vicinity of Tc for samples with the op- 
timal shape [Lt = L^"'^) to keep the second argument of 
the scaling function (0J constant. Fig. 01 shows a scaling 
plot ofgav vs. T for impurity concentration p = ^. Again, 
the data scale very well, but since the scaling function 
lacks the characteristic maximum, the error of the result- 
ing scaling factor is somewhat larger (1 . . . 2%) than 
that of Uf^"^^. The same quality of scaling was achieved 
for the other dilutions. Fig. shows the scaling factor xl 
vs. L for all four data sets. A combined fit to the ansatz 
XL — cL^/'^{\ -\- dL~'^^) where v and uj2 are universal, 
gives v = 1.16(3) and W2 = 0.5(1). As above, the fit is 
robust and of high quality (x^ « 1.2). Importantly, as 
expected for the true asymptotic exponent, v fulfills the 
Harris criterion 0i > 2/d=l. Note that both irrele- 
vant exponents loi and UJ2 agree within their error bars, 
suggesting that the same irrelevant operator controls the 
leading corrections to scaling for both z and ly. 

We have also calculated total magnetization and sus- 
ceptibility. The corresponding exponents P/v = 0.56(5) 
and 7/1/ = 2.15(10) have slightly larger error bars than z 
and I'. Nonetheless, they fulfill the hyperscaling relation 
2/3 + 7 = {d + z)v which is another argument for our 
results being asymptotic rather than effective exponents. 

In summary, we have performed Monte-Carlo simula- 
tions of a 3d classical Heisenberg model with linear im- 
purities which is in the same universality class as the 
dimer-diluted bilayer quantum Heisenberg antiferromag- 
net. We have shown that the QPT in this system is 
controlled by a conventional, finite-disorder critical point 
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FIG. 5: Scaling factor vs. L for four disorder concentrations 
p = i, i, f and i. Solid lines: Fit to xl = cL^/''(l + dL""^) 
with V = 1.16(3) and UJ2 = 0.5(1). 

with power-law dynamical scaling and universal expo- 
nents. (Note that the Ising version of our model, the 
diluted 2d random transverse I sing model, shows an 
infinite-randomness critical point |llLl28l| .l 

Let us compare our results to previous work. The mul- 
ticritical point &t p — pp and ~ 0.16, found in 

Refs. |ii|23, has a dynamical exponent z « 1.3. Within 
the error bars, this value coincides with the one found 
here for the generic p < Pp transition. We see no a-priori 
reason for this coincidence, so far it is unclear whether 
or not it is accidental. Vajk and Greven also quote 
exponents for p < pc- At dilution p = 0.25 they find 
z = 1.07 and v = 0.89, different from our results. How- 
ever, as the authors of Ref. ^3 point out, a value of 
V < 1 violates the Harris criterion, indicating that it 
represents an effective rather than an asymptotic expo- 
nent. It would also be useful to compare our findings 
with analytical results. To the best of our knowledge, 
the only quantitative result is a resummation of the 2- 
loop e-expansion |27l| . The predicted exponents signifi- 
cantly differ from ours; but they also violate the Harris 
criterion, casting doubt on their validity. 

Finally, we comment on experiments. If chemical dop- 
ing replaces magnetic by non-magnetic ions in an antifer- 
romagnet, e.g., Cu by Zn in YBa2Cu306, the case of site 
rather than dimer dilution is realized. The most promis- 
ing way to achieve bond dilution is the introduction of 
strong antiferromagnetic intra-dimer bonds at random 
locations. Thus we propose to study magnetic transi- 
tions in bond-disordered systems; those transitions can 
be expected to be in the same universality class as the 
one studied here. One candidate material - albeit 3d - 
is (Tl,K)CuCl3 [i^ under pressure; interesting quasi-2d 
compounds are SrCu2 (603)2 or BaCuSi206, where suit- 
able dopants remain to be found. 
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